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I . INTRODUCTION 

The period covered by this report has been spent primarily in coding and 
verifying the equations of motion for the Earth-Moon system as presented in 
Reference [1] . Some attention has also been given to numerical integration 
methods and parameter estimation methods. Existing analytical theories such as 
Brown's lunar theory as updated in [2]; Eckhardt's theory for lunar rotation, 

[3] ; and Newcomb's theory for the rotation of the Earth, [4], have been coded 
and verified. These analytical theories serve as checks for the numerical 
integration. Laser ranging data for the period January 1969 - December 1975 
has been collected and is stored on tape. 

This report presents descriptions and verifications of the several programs 
developed to date. A discussion of the parameter estimation method to be used 
and certain supporting theoretical developments are also presented. 

‘ : ■■ ■ ■ ' ■ ■ ■ ' ' I 

The main goal of this research is the development of software to enable 
physical parameters of the Earth-Moon system to be estimated making use of the 
extremely accurate data available from the Lunar Laser Ranging Experiment 
(LURE) and the Very Long Base Interferometry experiment (VLBI) of project 
Apollo. 

A more specific goal is to develop software for the estimation of certain 
physical parameters of the Moon such as the inertia ratios a , 8 , y , and 

the third and fourth harmonic gravity coefficients, C 1 3 j , S' 3 ^ , C'i^ , 

S '4 (j = 0, 1, 2, 3; k = 0, 1, 2, 3, 4), A unified model of the translational 

It 

and rotational motion of the Moon is to be utilized in the estimation process. 
Also LURE data only will be processed. 

In preparing the software, several basic questions must be asked, viz.‘, 

A. Should existing software be used to model the dynamics of the Moon? 

B. What type and order numerical integration scheme should be used? 

C. How are the numerical integration results to be checked? 

The answer to these questions — in part — is given below: 

A. Existing software provides a piecemeal treatment of the problem in that 
three separate numerical integrations must be made: 


i) Integration of all planets and Earth-Moon bary center, 

ii) Integration of Moon with respect to Earth-Moon barycenter. 

iii) Integration of lunar rotational equations. 

This is the approach used by the LURE Team and presented in the literature, 
[5] . This approach ignores the coupling of translational and rotational motions 
when terms of order (l/r 1 *) in the mutual distance are retained in the differential 
equations. It also does not provide for the inclusion of mutual gravitational 
potential terms which arise when terms of (1/r 5 ) are retained. Terms of at 
least the above orders must be retained to give physical libration accuracies 
of .005” to .01". This corresponds to a 2 cm to 3 cm accuracy in the LURE data» 

In order to achieve a unified numerical model and to allow rigorous 
investigations of the coupling and mutual potential effects it was decided to 
develop new software. The dynamical model was discussed in [ 1] and is referred 
to as a unified model of lunar translation and rotation (UMLTR) . it is a 
numerical integration of the combined lunar translational and rotational 
equations providing rigorously for coupling and mutual potential terms. This 
model is an integral part of two programs to be described later, viz., RIGEM 
and ESTEM. 

B. The success of Oesterwinter and Cohen [6] in integrating the solar 
system using a high order "Cowell" type method led to the development of a 
numerical integration subroutine based on that method. In the process of 
developing and verifying this routine, a numerical integration developed by 
Everhart, [7], RAI9S was obtained. This is equivalent to a high order 
implicit Runge-Kutta scheme and is an implicit single- sequence method using 
Gauss-Radau- and Gauss-Lobatto spacings. A new "Cowell" type routine started, by 
RAI9S was then developed and is referred to as COW. Subroutine RAI9S has been 
used in all runs to date in the development of the dynamics although a comparison 
■of RAI9S and COW is anticipated. The best integration order should also be 
decided based on a future study. 

C. The practice in the literature has been to compare numerical models 
of lunar rotation with existing analytical theories. Eckhardt’s theory [3] as 
modified by Williams [8] has been the basis for comparison. More recently, 
Eckhardt has developed a "300 Series" [9] lunar libration model. 
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To verify the unified model, comparisons must be made with both a lunar 
physical libration model and a solar system integration model. Initial 
conditions for the integration must be determined to obtain the best fit 
between the theory and the integration. 

Program ANEAMO was developed to provide an evaluation of existing lunar 
libration theories for comparison purposes. The solar system motion to date 
has been compared with data found in Oesterwinter and Cohen's work, [6]. 

In summary, the programs developed or utilized to date and their 
capabilities are: 

ANEAMO (01): This program evaluates 1) a truncated form of Brown's lunar 

theory, 2) Eckhardt's lunar physical libration theory, and 3) Newcomb's theory 
for the rotation of the Earth. It provides printed or punched output for use 
in other programs. 

RAI9S: A single sequence integrator developed by Everhart, [7]. It 

integrates a system of N first or second order equations with orders of 7, 

11, 15 or 19. Another similar routine RADAU31 has also been checked out which 
is capable of orders 7, 11, 15, 19, 23, 27, and 31. This latter routine 
requires double precision computations on the CDC 6600 series computers however. 

COW: A multistep predictor-corrector method for the solution of first or 

second order systems of differential equations, [101. The method is based on 
Bessel's central difference interpolation polynomials and the assumption that 
differences of some even order remain constant. The method is started with 
RAI9S. Orders of 4, 12, and 16 have been utilized. Starting is accomplished 
with a comparable order of RAI9S. 

RIGEM (01) : This program numerically integrates the combined translational 

and rotational motions of the Earth and Moon together with the translational 
motion of the remaining planets and Sun treated as particles. Options exist 
for integrating a subset of the planets and for eliminating Earth rotation. 
Subroutine RAI9S is utilized. Printed output or punched output is available. 

ESTEM: The basic putpose of this program is to fit 'the numerically 

calculated values of the lunar libration angles using the UMLTR to those obtained 
from an analytical theory. The fit is made using an iterative weighted least 
squares approach with a , 3 , y , and the initial Euler parameters and their 
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rates being adjusted in the process. The basic dynamic model is RIGEM (01) 

(with Earth rotation integration removed) . Necessary partial derivatives 
are generated using a numerical "variant trajectory" approach. The final 
residuals of the fit are output in plotted form. 

This report includes documentation for ANEAMO (01) , RIGEM (01) , and ESTEM 

( 01 ). 

Finally a verification run of ESTEM is described which compares the 
UMLTR with Eckhardt's libration theory for one year. Final residuals are 
periodic with amplitudes of ±3" in the latitude librations la and p and 
±10" in the longitude libration, t . 

Work is in progress to reduce the size of those residuals to the 1" level. 

Listings and decks for all programs and subroutines are available from the 
author. 

II. PROGRAM DESCRIPTION AND VERIFICATION 
A. Program ANEAMO (Version 01) 

General. Program ANEAMO (01) provides the capability for evaluating 
analytical theories of the translational and rotational motions of the Earth 
and Moon. Version 01 provides the capabilities for 

i) Evaluating a truncated form of Brown's lunar theory, 

ii) . Evaluating Eckhardt's theory for lunar physical librations. 

iii) Evaluating Newcomb's expressions for the precession, nutation, and 
spin of the Earth, 

iv) Providing orientation information in several useful forms such as 
direction cosines, Euler parameters, and the axis and angle of rotation, and 

v) Providing punched card output of lunar physical librations in 
longitude and latitude. 

Program and Subroutine Descripti on. Salient features of the main program 
and subroutines including theoretical developments are presented here. Further 
information may be found in the program listings in Appendix A. Refer to 
Reference [1] for general theory. 

9 

Program ANEAMO . This program handles all input/output operations. It 
calls in succession the subroutines necessary to evaluate the various theories. 
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The program also provides for the calculation— by differentiation of a cubic 
spline fit — of the time derivatives of the Euler parameters, (3 V) and {$'" } as 
evaluated from the rotational theories. These may then be used as initial 
guesses at the initial conditions for numerical integrations of the rotational 
equations of motion. The spline fit is accomplished using a standard subroutine 
SPLDER, [11] . 

This program evaluates all quantities at a series of Julian dates from some 
initial date VJIN to a final value VJF in steps of VJINC days. 

A multi-case option may be exercised if the parameter ICODE = 1, otherwise 
if ICODE = 0 then the program stops. 

Finally, the output parameter IPT may be set equal to 1 if it is desired 
to have the physical librations in longitude and latitude punched on cards. If „ 
this option is selected, a printout of these quantities is also made. 

Program ANEAMO links several subroutines as follows: 



Subroutine Luth . This subroutine contains a truncated form of Brown's 
lunar theory taken from Reference [2]. The fundamental arguments: 

(T , mean longitude of the Moon (VI) , 
r , Sun's mean longitude of perigee (V2) , 

P' , mean longitude of lunar perigee (V3) , 

ft , longitude of mean ascending node of lunar orbit on the ecliptic (V4) , 
D , mean elongation of Moon from Sun (V5) , 
e , mean obliquity of the ecliptic (V6) , 

were programmed as they appear in Reference [12] . They are equivalent to 
corresponding arguments L , O' , S , ft , D , e of the Improved Lunar 
Ephemeris, [2]. 
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The Moon's longitude, latitude, and parallax with respect to the mean 
equinox and ecliptic of data are calculated from the formulae: 


Longitude = L (i6) + 6 l (in) + <$L (ia) + 6L (i 6) (1) 

S = F (i0) + SL (in) - <Sft (in) + 5S (ig) (2) 

Latitude = A sin S + B sj.n 3S + C sin 5S + DN (ig) 

+ 6 Lat (ie) (3) 

sin II = (5 sin IT (iy) + 5 sin H (i?) ) (1 - 4.6747 * 10“ 5 ) (4) 

Parallax = sin II + (1/6 sin 3 (II) )/ (206265) 2 (5) 


where 

L (i0) = <r 

F (i0) * (T - Si 

SL (in) = additive terms in longitude 
<$L (ia) = solar terms in longitude 
6L ( i5 ) = planetary terms in longitude 
(Sfl (in) = additive terms in node 
6S (ig) = solar terms in latitude 
5 Lat (ie) = planetary terms in latitude 
N (ig) = solar terms in latitude 

a = yi + Y] c ■ ;; . . 

B = Y 2 /Yl A 
c = Y 3 /Y 1 A 
D = 1/Yl A 

Yi C = solar terms in latitude 
Yj , Y 2 • Y 3 = coefficients of principal terms in latitude. 

The principal terms in latitude Yl » Y 2 * Y3 were respectively multiplied 

by the first, third and fifth power of the factor (1 + 2.708 x 10" 6 + 139.978 Sy ) 

C 

where Sy^ are additive terms from list in *■ 
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All notations used above are explained in the ILE. Certain small 
correction factors have been ignored that are not consistent with the accuracy 
of the terms retained. Table 1 presents the accuracy of terms retained in 
version 01 of ANEAMO. 

Table 1. Truncation Level of Brown's Lunar Theory. 


Coordinate 


Amplitude of smallest terms 
retained (rounded to value 
shown) 


Longitude 

Latitude 

Parallax 


Additive terms due to planetary perturbations and other small long period 
effects are included in (£ , ft , and 6y c . Table 2 lists the terms programmed 
in version 01 of ANEAMO. 


Table 2. Additive Terms Programmed in ANEAMO. 


Ref. No. 

Add To 

Coefficient 

1365 

L 

0.84 

1366 

L 

0.31 

1373 

L 

14.27 

1375 

L 

7.261 

1376 

L 

0.282 

1379 

L 

0.075 

1382 

L 

0.237 

1383 

L 

0.108 

1385 

L 

0.126 

1369 

ft 

0.63 

1406 

ft 

0.17 

1407 

ft 

95.96 

1408 

ft 

15.58 

1409 

ft 

1.86 

1413 

6 V 

-4.318 

1414 


-0.698 

1415 

6y c 

-0.083 


Planet 


Venus 

Venus 


Mercury 

Venus 

Venus 

Venus 


Venus 















X 


I 
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This subroutine also calculates the geocentric ecliptive coordinates 

{x.} and geocentric equatorial coordinates (x. } of the Moon. These are 
X ECju X EQ 

mean equinox and ecliptic of date and mean equator and equinox of date systems 
respectively. Formulae utilized for this computation are 


r = (6378. 16) /Parallax 


* X i*ECL 


r cos (latitude) cos (longitude) 
r cos (latitude) sin (longitude) 
r sin (latitude) 


{x i } EQ 


'1 0 

0 cos e 

.0 sin e 


0 

-sin e 
cos e « 


{x t } ECL 


(6) 

(7) 

(8) 


Subroutine LURO. This subroutine calculates the physical librations in 
longitude, node, and inclination based on Eckhardt's librational theory contained 
in References [3] and [13] . The additive and planetary terms constructed by 
Williams in Reference [8] are also included. This is basically a second degree 
theory in the lunar gravity harmonic coefficients. Eckhardt's theory provides 
the physical librations in the form 


where 


T. = V T. . sin (NL. • & + NLP. • V + NF. • F + ND. ® D) 

3 fr* XD x x x x 

x 


Id. = YL C t, • sin (NL. * 1.+ NLP « + NF • F + ND » D) 

3 ^ K3 Jc jc jc Jc 


(9) 


p . = V R . sxn (NL • £ + NLP 0 SL + NF • F + ND « D) 
3 mi m m m m 

Kl 



physical libration in longitude 
physical libration in latitude 
physical libration in inclination 
inclination of lunar equator to ecliptic 
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j = index number for set of coefficients being used 


i, k, m = summation indices 
l = L - Si 
£' = L' - S' 

f = l - a 

D = L - L' 

NL, NLP, NF, ND = integer coefficients 


The coefficients 
of the Moon, where 




R . depend on the inertia ratios 0 
W3 


Y 


0 = (C - A)/B 
Y = (B - A)/C 


and A , B , C are the lunar principle inertias. Thus, the j index 
distinguishes between the coefficients for different values of 0 and y . 
Currently coefficients for three j's have been used as shown in Table 3 (see 
Appendix A for numerical values) „ 


Table 3. Eckhardt Coefficients Programmed. 


j 

I 

6 

Y 

1 

5521.6" 

6.268 x io“ 4 

2.3 x lO -4 

2 

5559.6" 

6.3 x 10 -4 

2.0 x lo” 4 

3 

5550.2" 

6.293 x 1Q~ 4 

2.27 x i(T 4 


The additive and planetary contributions as calculated by Williams in [8] 
are for 0 = 6.293 x 10~ 4 and y * 2.27 x 10 -4 . The j = 1 coefficients 
come directly from one of Eckhardt 's tables and are rounded to 1", The 
coefficients in j = 2 and j = 3 were obtained by interpolating between various 
tables given by Eckhardt keeping terms to 0.2" in one of the tables interpolated 




The inclination of the lunar equator to the ecliptic, 1 , is also a 
function of 3 and y . The values used for I are also shown in Table 3 
as calculated from Eckhardt ' s empirical formula 

I = -1612“ - 5.2 x 10 4 y + 11.4 x 10 6 8 . ( 

Once the physical librations are calculated, the Moon's orientation in 
inertial space is calculated as follows. 

The physical librations T , Ic , p . are perturbations to three Euler 
angles ip , 0 , <f> locating the Moon with respect to the ecliptic and mean 
equinox of date. To express this mathematically, consider the geometry of 
Figure 1. The {z^} axes are fixed to the rigid Moon and are principle center 
of-mass axes. The axes {x. } are the ecliptic and mean equinox of date 
set. The axis notation used here follows that defined in Reference [1] . 



zi 


Figure 1. Orientation of Moon in Space. 



/ 


The above axes are related as follows; 


(z.) = [KM] {x i ) EcL 


" (o/i op - Sp S\p C0) (C4> Si|/ + S<f> C0 Op) -Sp S6 

[RM] = (-Cip Sp “ Cp C0 Sp) (-Sp Sp + C? Op C0) -Cp S0 

-Sip S0 Op S0 C0 


The matrix [RM] is a rotation matrix for a rotation from {x. , to {z.} 

i ECL i 


in the sense ZXZ through angles ip , -0 , 

S<}> E sin “p and Cp = cos 'p is used above. 


The shorthand notation 


Now, the Euler angles and physical librations are related through the 


expressions 


0 = I + p 
ip = ft + a 
$ = h + (T-iP + t 


where I = constant and ft and are the fundamental arguments from the lunar 
theory. j 

The orientation of the Moon with respect to the mean equator and equinox of 
1950.0 frame, {X^ 1 } , may be obtained as follows: 

{zj = [RM] [EC] T [P] {X i 1 } - (13) 


where 


{x.} eq - [EC] {*.} ECL 

{x.) eb = m {x. •> 


[P] = Precession matrix (defined later) 
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. ’1 0 o c 

[EC] = 0 ce -se , 

.0 SG CE . 


The orientation of the Moon with respect to the reference axes {Z^ for 
the Moon may be obtained as follows: 

{zj = {XM<|>] {Z ± } . (14) 

In order to evaluate matrix [XM<{>] the relation between {X^' } and {Z^} 
must be known. This depends only on the position of the Moon's center of mass 
with respect to {X^} as discussed in Figure 2 of Reference [1], Referring to 
that figure, the relation between {X^' } and {Z^} can be obtained by 

i) Rotating {x. 1 } through (X + II) about X 3 ' to obtain {Z.'} 

JL 1 

and then 

ii) Rotating {z^‘ } through 4> about Z 2 ' to obtain {Z^} . 

The compound rotation provides 


{Z i > = [T] {X i '} 


where 


[T] = 


■ -CX C<J> 

-C<p sx 

-S<J) 

sx 

-CX 

0 

. -CX s<p 

-SX S(j) 

C<t> 


Here 4> and X are the geocentric latitude and longitude of the Moon with 
respect to the mean equator and equinox of 1950.0. The Moon's position, ^ X i^EQ 50 
is obtained from Brown's lunar theory as calculated in LUTH. Thus, 


{x i } EQ50 “ IP1 {x i } EQ ’ 


Next, the polar coordinates r , , X as required in (15) are found from 
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r = /x? 


+ x 


+ xi 


C<J» CX = xi 


1EQ50 2EQ50 3EQ50 


EQ50 / 


C * CX = X 2 E Q50 /r 
s <t> = x 3 EQ50 / r 


C * = /x lE 2 50 + X 2EQ50 /r 
SX = C4> SX/C<J> 

CX = C< J> cx/c<t> . 


( 16 ) 


Finally, the matrix XM(j> is constructed as follows 


[XM<f>] = [RM] [EC1 T [P] [T] T , 


(17) 


thus providing the orientation of the Moon with respect to the {Z^} reference 
axes . 

LURO also calculates the Earth's selenographic coordinates, X_ and . 

By definition, the Earth's selenographic coordinates are the latitude (y e ) 
and longitude (X ) of the Earth as seen from the { z .} axes. Thus, 

© X 


Cp e cx e ~ ~ k l " i r = + *1 • Ki 


cu e sx e = -k 2 a i = k 2 * K 2 


(18) 


sy e = -k 3 • i r = k 3 * k 3 


where i is a unit vector directed from the Earth' s mass center to the lunar 
r ' 

mass center as shown in Figure 2 of Reference [1] . Since i^ = -Ki , the 
direction cosines in (18) are just the elements of the first column of matrix 
XM<f> . Then 


X = tan' 


_1 XMtfr (2,1) 

XM$ (1,1) 


Vi = tan 


-1 


XM4> (3,1) 


/ XM(J) (1,1) 2 + XM<J> (2,1)2 J 


(19) 


13 


) 


Subroutine EARRO . Version 01 of EARRO contains the approximate expressions 
for Earth rotation provided in Reference [4] . 

The true sidereal time, accurate to 0.2" or 10“ 6 is given by 

0 =* 100.075542 + 360,985647348 T 

+ 0.29 x 10~ 12 T 2 - 4.392 x 10~ 3 sin (ax) 

+ 0.053 x io~ 3 sin (a 2 ) “*0.325 x io~ 3 sin a 3 

- 0.05 x io“ 3 sin (a 4 ) ' (20) 


where 


T = Julian Date - 2433282.5 
a! = 12'? 1128 - 0,052954 T 

a 2 = 2a i 

a 3 = 2 (280?0812 + 0.985647 T) 
a 4 = 2 (64?3824 + 13.176398 T) . 

The transformation from the {x^} system to the {y^ system is given in 
Reference [4] as 


{y ± } = [S] [N] [P] {x i 1 } 


(21) 


where [S] -is the spin matrix, [N] is the nutation matrix, and [P] is the 
precession matrix. 

More explicitly. 


[S3 


[P] 


’ ce se o' 

I -se ce o , 

.o o i . ■ ; 

' (-Sk Sco + Ctc Cto CV) 
(Sic Cco + Ctc Sco Cv) 
CK SV 


(“Ctc SO) - SK Qi) CV) I -Cut SV 
(Ck Cco “ Sk Sco cv) j -Sco Sv 
-Sk Sv | cv 


( 22 ) 


(23) 
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and 


[N] = 


C (Av!) C (-Ay) I 

(C (-Am) s (-Ae) S (Av) 1 
- s (-Am) C (-Ae)) ! 

(C (-Am) S (Av) C (-Ae) | 
+ S (-Ae) S (-Am)) I 


C (AV) S (-Am) I - s (AV) 

(S (-Am) S (-Ae) S (Av) 1 S (-Ae) 

+ C (-Ae) C (-Am)) | C (Av) 

(S (-Am) S (AV) C (-Ae) j C (-Ae) 

- S (-Ae) C (-Am)) I C (Av) 


(24) 


The rigorous form of the nutation matrix rather than the approximate form 


IN] - 


1 

am 
L Av 


-Am 

i 

Ae 


-Av 

-Ae 

1 


(25) 


was programmed to insure that [N] was rigorously orthonormal. 
The arguments utilized in the above matrices are 

* 

k = 0S063107 T 
co = 0?063107 T 
V = 0?0548757 T 

Am = -76.7 x 10 -6 sin (a^) + 0.9 x 10“ 6 sin (a 2 ) 

-5.7 x io“ 6 sin (33) - 0.9 x io~ 6 sin (a^) 

AV = -33.3 x 10~ 6 sin (aj) + 0.4 x 10" 6 sin (a 2 ) 

-2.5 x io~ 6 sin (33) - 0.4 x 10“ 6 sin ( a 4) 


Ae = 44.7 x 10“ 6 cos (a2) - 0.4 x 10“‘ b cos (a 2 ) 
+ 2.7 x 10~ 6 cos (33) + 0.4 x 10“ 6 cos (34) 


-6 


(26) 


(27) 


Note that Am , Av , Ae are given in radians. 

Reference [4] lists the accuracy of the above expressions as 0.2" or 
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-6 


The values given in equations (26) and (27) are programmed. If more accurate 
values are required, the following may be used. 
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K = (2 3042V 5 3 + 139V73 T + 0V06 T 2 ) t 
+ (30V23 - 0V27 t) t 2 + 18V00 t 3 

u - K + 179" 27 + 0V66 t) t 2 + 0V32 t 3 (28> 

V = (20046V85 - 85V 33 T - 0V37 T 2 ) t 
+ (-42V67 - 0737 t) t 2 - 41.80 t 3 


where r is the epoch (1950.0) of {X. ' } and t is the epoch of the mean 

i 

sidereal system, both measured in thousands of tropical years from 1900.0. 

Also many additional nutation terms are listed in References [2] and [14] . 

Subroutine ICOND . This subroutine accomplishes a transformation from the 
Euler parameters { 8 "} for the Earth to the Euler parameters { 3 ' } . The 
relation between those was presented in equation (16) of Reference [1] and is 

{ 8 ' } = 18 ] ^ { 8 ”} - 


The Euler parameters ( 8 ) locate the reference axis system {Y^} with 
respect to the system {x^ 1 } . Since the {Y^} system is rotated, through an 
angle, a , about the 13 ' axis with respect to (x^ 1 } can be shown that 

80 = cos (a/ 2 ) 

81 - 0 

(29) 

82 - 0 

83 = sin (a/ 2 ) , 

where 


a = ao + a T. 

Thus, 


I 


c (a/ 2 ) 

0 

0 

-S (a/2) 

: j 1 

0 

C (a/2) 

-s (a/ 2 ) 

0 

j 

■ ; i . 

0 

S (a/2) 

C (a/2) 

0 


s (a/ 2 ) 

0 

0 

c (a/ 2 ) 

• 



(30) 


( 31 ) 


Due to the nature of [8] 

ter 1 - [0] T . (32) 

Currently , 

a Q = 100.075542 
a = 360.985647348 ; 

are programmed. 

Subroutine AXANG . This subroutine computes the axis and angle of rotation 
from any rotation matrix [R] as well the corresponding Euler parameters 
{6^} . The sense of the rotation is from (x) to {x' } where 

{X*} = [R] {x} . (33) 

The formulae for this are provided in Reference [15] and are summarized below. 
The angle of rotation is 5 and the direction cosines of the axis--with respect 
to both {x} and {x' } — are {c^} - 

cos 5 = 1/2 (Rji + &22 ^3 3 ~ 1) 

C 1 = ( a 23 “ a 32)/ 2 si 11 5 

c 2 “ ( a 3 i “ aj 3)/2 sin 6 ■ 

C3 = (a 12 “ a 2 1 ) / 2 sin 5 • 

The Euler parameters for the rotation may be found from 

So = cos (6/2) 

S. = C. sin (6/2) 

1 1 

logic is incorporated in this subroutine that keeps the calculated rotation 
axis generally aligned with the body rotation axis. 


(36) 


(i = 1, 2, 3) 


(34) 

(35) 
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Verification of ANEAMO. The verification of program ANEAMO has been 
accomplished in several ways 

i) Table 4 provides a comparison of the calculated values of the fundamental 
arguments of the lunar theory as calculated by ANEAMO and as given in the American 
Ephemeris and Nautical Almanac (AENA) for 1974, [12} . 

ii) Table 5 provides the residuals in geocentric ecliptic longitude, 
latitude and parallax from a comparison of ANEAMO with the AENA. The sense of 
the residuals is ANEAMO-AENA. Note that the nutation in longitude must be added 
to the values calculated in ANEAMO for the comparison. The calculations are made 
for a two-month period beginning at J. D. 2442050.5. Table 5 shows actual residuals 
residuals, their mean (x) and standard deviation (5D) for the case when the 
longitude and latitude series are truncated at 1" and the parallax series is 
truncated at O’. '01 . Additional terms in the series axe to be programmed but the 
accuracy shown in Table 5 is sufficient for present purposes. 

iii) The sidereal time, 0 , as calculated in EARRO was compared with the 
value given in the AENA for J. D. 2442050.5. A residual of 0'.'3 resulted. 

iv) The calculation of the nutation matrix elements Av , -Ae , and 
-Ay was compared with values presented in the AENA. For J. D. 2442050.5, 

ANEAMO provides. 

Ay = 0.78227 x 10~ 4 rad. 

Av = 0.33966 x 10“ 4 rad. 

Ae = -0.36289 x 10“ 5 rad. = -0.748" 

Comparable values to the above are not presented in the AENA but Aip and Ae 
are given there. The values are 

Ae = -0V740 

A^ = 17'.'493 . 

The nutations Ay and Av are related to Aijj through 

Ai|) = Ay/cos e = Av/sin e 
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Table 4. Verification of Fundamental Argument 
Calculation in ANEAMO 


Quantity 


tt 

++ 

ANEAMO 

AENA 

35? 02006 

35f020l 

282?49337 

282 ?49335* 

105? 63199 

105?6320 

267 ?81342 

267?8134 

112?78291 

112?7829 

23.°44266 

23?44266 + 

289?38807 

289?3881 

282?23714 

282 .“2371 

359 “74377 

359?74375 

127.°20664 

127.“2067 


All comparisons made for J. D. 2442050.5. 


* Hand calculation. 


Tabular value of e in AENA is an "of date 
value. This is related to the mean value, 
T , as calculated in ANEAMO by e 0 D ® ” 

+ Nutation in obliquity. u 




















thus the residuals in Ae and Ai|> may be calculated as 

Res. in At|> = 0'.'1 
Res. in Ae = 0'.'008 . 

v) The precession matrix [P] as calculated .in ANEAMO for the date 1974.0 
(J. D. 2442048.2358) could be compared with the same matrix appearing in the 
AENA, Out of the nine elements the worst case residual was 4.76 x 10“ 7 . 

vi) The verification of the calculation of p , la , and T , the 
physical librations, is manifested by a comparison of the theoretical values 
computed in ANEAMO with numerically integrated values. This comparison will 
be discussed later. 

vii) Any set of Euler parameters must satisfy the following constraints: 

S 0 2 + 6i 2 + 3 2 2 + B 3 2 - 1 (37) 

3o 3 0 + 3i 3x + 3 2 3 2 + 33 33 = 0 (38) 

These constraints are tested in program ANEAMO. The constraint (37) is satisfied 

more accurately than (38) since in the latter the rates 8. are calculated from 

■ 1 

a cubic spline fit to the parameters 8^ . Typical values encountered are 
Moon: 

Z 6. 2 - 1 

2 3. 3. 

1 1 

Earth : 

2 3. 2 - 1 = 3.2 x 10” 14 
x 

2 3. 3. = 2.2 x icr 8 
xx 

B. Subroutines RAI9S, RADAU31 

These are numerical integration routines for systems of first or second order 
ordinary differential equations. The theoretical development of the method used 


= 1.8 x 10~ 14 
= -7.5 x icr 13 
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is presented in Reference [7]. Basically the solutions to x.= F (x, t) are 
developed in truncated series in time t whose coefficients are found empirically. 

The method is a single-sequence method that uses Gauss-Radau and Gauss-Lobatto 
spacings for the several substeps within each sequence. The method is equivalent 
in principle with the implicit Runge-Kutta methods. 

Subroutine RAI9S . This is a single precision deck suitable for a computer 
with a 60 bit word length in that precision. Integration orders of 7, 11, 15, 
and 19 are provided. This routine is used both in programs RIGEM and ESTEM and 
as a starter to the Cowell second-sun method of subroutine COW. 

Subroutine RADAU31 . This is a double precision deck suitable for CDC 
machines with a 120 bit word or for IBM machines with a 128 bit word. Integration 
orders of 7, 11, 15, 19, 23, 27, and 31 are provided. 

Verification. Before use both of the above subroutines were verified by the 
calculation of a test orbit from the restricted three body problem, [10] . The 
orbit is shown in Figure 2. It is a periodic orbit in the rotating frame with 
period 

t f = 6.1921693313 19639 70699 23217 ... 

Initial conditions are 

Yl = 1.2 yi' » 0.0 

y 2 = 0.0 

y 2 ’ = -1.04935 75098 3031990731 0410434 . „ . 

The orbit has three loops and requires frequent step size changes. 

Table 6 presents the results of the verification. 

C. Subroutine COW . 

The success of Oesterwinter and Cohen in integrating large systems of 
equations using high-order multistep methods [6] led to the development of this 

■ • ■ . ■ ' . ; j \ ' j ... 

subroutine. The method was programmed in accordance with Reference [10], It is j 

referred to as Cowell's method and is a multistep predictor- corrector routine for 
systems of first or second order ordinary differential equations. Cowell's method 
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Figure 2. Test Orbit for Subroutines RA19S and RADAU31 
from Restricted 3-Body Problem. (^ = Earth 
(£ = Moon) 


Table 6. Accuracy of RAI9S and RADAU31 on test orbit. 


Method 

Order 

Ayi 

! 

Ay x 

Function Calls 

RAI9S 

7 (5)* 

1.7 x 10^ S 

3.4 x kT 5 

736 

RAI9S 

11 (6) 

1.1 x io~ 7 

1.9 x lo~ 7 

1025 

RAI9S 

11 (7) 

4.6 x 10~ 9 

7.0 x lo” 9 

1437 

RAI9S 

15 (10) j 

4.5 x 10“ 13 

4.9 x 10~ 13 

2867 

RAI9S 

19 (12) 

2.3 x 10“ 15 

2.6 x 10“ 15 

3802 

RADAU31 

23 (15) 

• 2.1 x 10” 21 

3.6 x 10“ 19 

5820 

RAD AU 31 

27 (20) 

1.0 x 10~ 24 

1.6 x 10~ 24 

39610874 


* Numbers in parentheses are sequence sise control numbers. 
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is based on Bessel's central difference polynomials and the assumption that 
differences of some even order iremain constant, [10] . The coefficients fof 
4th, 12th, and 16th order methods have been punched. The method is started 
using RAI9S ojf a comparable order. A constant or variable step capability as 
well as an "exact end" capability is included. The "exact end" capability uses 
the RAI9S Routine also. Options exist for no corrector, n applications of 
the corrector or iteration using the corrector until convergence is obtained 
within a specified accuracy level. 

Verification. The same test orbit utilized for RAI9S has been integrated 
using COW. The results are shown in Table 7. 


Table 7. Verification of COW using test orbit. 


■ 

Order 

Halving/Doubling 

Limit 

Ayi 

Ayi 

Function Calls 

Start 

Run 

11 

12 

1. x 10~ 10 

2.0 x 10“ 7 

3.3 X ICT 10 

8265* 

11 

12 

i. x ict 11 

2.2 x lo“ 10 

1.4 x 10“ 10 

10180 

11 

12 

5. x icr 12 

3.8 x 10“ 9 

1.2 x IQ- 8 

10761 


* Iteration using corrector formulae produced large number of function calls. 


The orbit integrated here places a stringent test on COW since step size 
changes are clumsy to handle with these types of methods . More work needs to be 
done on this subroutine in terms of optimization and criterion for choice of 
order and other parameters. 

D. Program RIGEM (Ve rs ion 01) 

This program provides the capability for numerically integrating the 
coupled translational/rotational motions of the Earth and Moon treated as 
arbitrary rigid bodies together with the remaining planets and Sun modeled as 
particles. The general theory is outlined in Reference [1] . Subroutine RAI9S 
is utilized as the integration routine. Options exist for 1) omitting Earth 
rotation; 2) omitting planets Mercury, Saturn, Uranus, Neptune and Pluto 
(<J>MPL (1) = 1); 3) for multiple cases; and 4) for punched output of the 
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calculated values of the lunar physical librations. The main program links 
several subroutines as follows 



Program RIGEM. This program calls subroutine PROB to obtain all initial 
conditions and values of parameters. It then provides an integration loop that 
calls subroutine RAI9S at TINC intervals from VJDEP to VJDEP + TMAX. The loop 
also provides for output of calculated quantities and for calculation of the 
Earth’s selenographic coordinates, the lunar physical librations and Euler 
parameter tests for Earth and Moon orientation. 

At each time step subroutine RAI9S returns the following current values 
to RIGEM: 

A. Current time TWR - VJDEP + TF . 

B. Position and velocity of planets and Sun with respect to the {x^' } frame. 

C. Euler parameters and rates {g^'} , { 0^ ' } , {g^ ,M } , {g^'" } ; as 
listed in Table 8. 

This program calculates the Euler parameter tests for the Earth mentioned 
earlier 

Z 8.‘ - 1 

; i 

z 6 i' : V “ 0 * , 

Since the immediate use of this program is in analyzing lunar motion, the 
lunar rotation segment has been more thoroughly treated than that of the Earth. 
Checks for gross errors have been made in the Earth rotation logic but no 
definitive verification studies have been made to date. Currently the Earth 
orientation is provided by subroutine EARRO as discussed earlier. 
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Translation 


Planet 

Mass 

Position 

Velocity 

Sun 

1 

x (1 - 3) 

v (1 - 3) 

Mercury 

2 

x (4 - 6) 

v (4 - 6) 

Venus 

3 

x (7 - 5) 

v (7 “ 9) 

Earth 

4 

x (10 - 12) 

v (10 - 12) 

Moon 

5 

x (13 - 15) 

v (13 - 15) 

Mars 

6 

x (16 - 18) 

| v (16 - 18) 

Jupiter 

7 

x (19 - 21) 

v (19 - 21) 

Saturn 

8 

x (22 - 24) 

v (22 - 24) 

Uranus 

9 

x (25 - 27) 

v (25 - 27) 

Neptune 

10 

x (28 - 30) 

v (28 - 30) 

Pluto 

11 

x (31 - 33) 

v (31 - 33) 


Rotation 



Parameter 

Rates 

Earth 

{ 0 i * > , {6V} 

x (34 - 37) 

j 

v (34 - 37) 


Noon 


{e^"} {ev"} 


v (38 - 41) 


x (38 - 41) 
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The rotation matrix [C (3'")] defined by 

{z i > = [C (3"' ) J {Z ± > (39) 

can; be evaluated at each step using {3^'"} . The matrix [C (3'" )] is of the 

form given in equation (72) of Reference [1]. As shown earlier the Earth's 

selenographic longitude and latitude, X and y , can be calculated from the 

6 6 

elements of the first column of [C (S'")] as done in equation (15). 

Next, the location of the Moon with respect to the Earth componentiated in 
the {X^ 1 } frame can be found as follows; 


DAL (1) = x (13) - x (10) ' 

DAL (2) = x (14) - x (11) ■ 

DAL (3) = x (15) - x (12) , 


(40) 


The polar coordinates r , X <j> may be found from equations (16) substituting 
DAL (1) for x 1 E q5q • etc * Next, matrix [T] of equation (15) may be formed. 

The physical librations, equation (12), may now be calculated as follows. 

The relation between { z ana {Z^} was derived earlier, viz. 


{ z ± } = [RM] [EC) T [P] [T] T {Z ± > . 


(41) 


Comparing equation (39) with equation (41) provides the result 
[RM] [ECJ T [P] [T] T = [C (3"’)] 
or 

(RM] = (C (3'")] [T] [P] T [EC] . : (42) 

Now [C (3'")] is available from RAI9S , [T] was calculated above, and [P] 
and [EC] are calculated as they are in EARRO. 

T 

The matrix [PP] defined by [PP] E [RM] is calculated in RIGEM, viz. 
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[PP] = [EC] T [P] [T] T [C (S'")) T . (43) 

The elements of [PP] (or [RM] ) are functions of 8 , if) , ip , the 
Euler angles locating the motion with respect to the ecliptic and mean equinox 
of date system. 

Accordingly , 


t = tan' 1 (-PP (3,1) /-PP (3,2)) 

0 = tan -1 ( / PP (3,1) 2 + PP (3,2) 2/PP (3,3)) 
ip = tan" 1 (-PP (1, 3)/PP (2,3)) 


(44) 


Knowing these angles the physical librations are 

p = 6 - I 

a = ip - ft ► 


(45) 


Finally, the Euler parameter constraints for the Moon are calculated, viz. 


2 3 . 


2 3 . 3 . 

i i 


2 3 . 3 . 

i i 


2 3. 2 


(46) 


Subroutine PROB, This subroutine provides all initial conditions and 
constants for RIGEM. These values are also printed out by PROB for reference. 
Initial conditions for the translational motion and mass parameters are currently 
being taken from Table 10 of Reference [6] . Initial Euler parameters and rates 
are taken from ANEAMO. 


If the multi-case option is selected, then a branch to statement 12 occurs in 
this subroutine. The initial values of the time, integration variables and 
their rates are reset automatically and whatever changes are desired in the 
parameters must be read in at this point. 
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Subroutine FORCE. This subroutine calculates the accelerations required by 
subroutine RAI9S. The calculations are made in the sequence shown in Table 9. 


Table 9* Calculation Sequence in Subroutine FORCE. 


Calculations 

Equations 
in Reference [1] 

A. F (4) - F (33) N Particle Accelerations on all 

10 

planets and Moon. 


B. F (1) - F (3) Accelerations on Sun. 

= 0 in Version 01 

C. F (34) - F (37) Euler parameter accelerations 

31 

{ fiT * } of Earth . 


D. F (38) - F (41) Euler parameter accelerations 

50 

{(3. } of Moon. 

. ■■ ■■ - --- „ , 



Basically (Version 01) subroutine FORCE calculates 

i) N Body gravitational forces on all bodies treated as particles based 
on equations (10) and (12) of Reference [1]. 

ii) Torques on the triaxial Earth due to a point mass Moon are given in 
equation (98) of Reference [1] . Note: Torques due to a point mass Sun are not 

included in this version. 

iii) Torques on the triaxial Moon due to a point mass Earth and a point 
mass Sun given by equations (84) and (91) of Reference [1] . 

All geometrical equations and transformation equations are given in equations 
(15) - (55) of Reference [1]; except those providing the location of the Sun with 
respect to frame {z^} . They will be developed in the next paragraphs. 

The torque components on the Moon due to a point - mass Sun are 
Mz 1 - 3 GM q a ihq n Q /r 3 0(L 

Mz 2 = -3GMq 3 *0 n Q /r 3 Q t (47) 

M 23 = 3GMq y *© n» G /r 3 ©^ 
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where , m^ , are the direction cosines of the Sun with respect to the 

{z.} frame. RIGEM provides the direction cosines of the Moon with respect to 
the {X.*} frame, viz. 



x (12 + i) 



(48) 


The direction cosines of the Sun with respect to a frame {x/ } at the Moon are 
the negative of the ratios given in equations (48) . 


The direction cosines of the Moon with respect to {z^} can now 
formed since the rotation matrices [T] and [C (8'")] are available (see 
equations (39) and (15)): 


{zA = [C ( 0 ** * ) ] izA 

{zA = [T] {X jL ’} . 

Accordingly, 

f~-] = -tC (S'" ) ] [T] I— ^ — } (49) 

Subroutine FORCE also calls subroutine FORTOR to provide the remaining 
forces and torques that are modeled. 

Subroutine FORTOR. This subroutine calculates: 

i) Force on Earth other than N particle force, 

ii) Force on Moon other than N particle force, 

iii) Torque on Moon other than gravity-: gradient effect due to Sun and 

Earth. 

iv) Torque on Earth other than gravity-gradient effect due to Moon. 

Specifically, in version 01, the following forces and torques are 
calculated : 
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i) Force on Earth — no additional contribution programmed, 

ii) Force on Moon — Force due to Earth's figure and Moon's figure, 

iii) Torque exerted on Moon because of lunar higher degree gravity 
harmonic coefficients and mutual potential terms by an oblate Earth, viz. 

a) Torque on Moon due to point mass Earth acting on Cgg' , C 31 ' » 

c 32' ' s 3l' ' £ 32 ' / c 3 3 ' ' s 33* ' c 4l' * S 41 ' , € 43 ’ , 
S 42 ' , and 

b) Torque on Moon due to the interaction of second degree lunar 
harmonics with Earth oblateness C 2 q C 2 o' > C 2 q C 2 2 ' • 

iv) Torque on Earth— -no additional contribution programmed. 


The above torque expressions were given in equations (85) and (90) of 
Reference [1]. 

The force expressions can be obtained from equations (59) - (78) of 
Reference [1] . They may be put in the following form: 


F = ~3G m 4 m 5 a" 
r 


r 45 

-3G HI 4 11150 a 12 
* 45 4 


[P20 C20 + ^21 (^21 CX + S21 SX) + P 22 (C22 C2 X + S22 S2 X) ] 
[C 2 0' p 20 + p 22 C 2 2' C2 X] 


— F, - — — [C 2 0 ^20^ + P 21x (C 21 CX + S 21 SX) + P 22 i (*-22 C2 X + S 22 S2 X) ] 

r v r 45 4 v 9 <p 

— ?4 — [O20 1 P20* + *22* C22’ 02 X] 

rits ■ 

X G m in 

F = [p 21 (_ C21 sX + S 2 i CX) + 2P 2 2 (-C 22 S2 X + S 2 2 C2 X)] 

rC<J> 


r45 4 C<j> 


,12 


1 G — [-2P 22 C 22 ' S2 X) 


r 45 4 C(j) 


(50) 


where r 4 5 , X , <j> are the polar coordinates of the lunar mass center with 

T 

respect to {z ± } . The quantities P 20 > P 21 r P 22 • p 20 ( j ) * p 2 li * p 22 

were defined in Reference [1) . 


* 
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The primed quantities refer to the Moon and the unprimed quantities refer 
to the Earth. Here, 


a 2 

M 

c 20 

= 1/2 [A + 

B + 

C - 3 

(a" 

' 2 A + 

S" 2 

B + Y" 2 C)] 

a 2 

M 

C 21 

= aa" A + 

S3" 

B + YY" 

C 




a 2 

M 

S 21 

= a' a" A 

+ 3* 

3" B + 

Y' 

y" c 



4a 2 

M 

c 22 

= A (a 12 - 

a 2 ) 

+ B (3 

12 

- S 2 ) 

+ c 

(Y 12 - Y 2 / 

2a 2 

M 

s 22 

= aa' A + 

33' 

B + YY' 

C 





where the relative orientation of the Earth and the Moon is given by the [A] 
matrix defined by 


{y i > = m {z ± } 

where 


(52) 


m 


a a' a" 

3 3 ' 3 " 

Y y * y** 


Now, EARRO provides (equation (21) ) 

{y ± } = [S] [N] [P] {X i ’} 

and RIGEM provides (equations (39) and (15)) 

(z. } = [C ( 3 '" ) ] [T] {X. * } 
i x 

so that 

[A] = [SI [N] [P] [T] T [C (S ,,, )] T . (53) 

The torques are programmed directly from equations (85), (89), (90). Note 
that the following simplifications have been made in the above torques as 
programmed in Version 01s 
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i) In the mutual potential terms, only those multiplied by C 2 o for the 
Earth are retained, 

ii) In the fourth degree terms only those factored by Cm* , 043 ' , 842 ' • 

and S 44 ' have been retained since these produce the most significant effects 
in the librations, [16]. 

The direction cosines of the Earth with respect to the {z^} frame as 
required by equations (89) and (90) are immediately available from the [C (3"')} 

A A 

matrix in RIGEM. The polar coordinates (r , <fi , A) of the lunar mass center with 
respect to {y^} are required in the mutual potential terms. These may be 
found as follows. 

Equation (21) provides 

{y ± > « [S] [N] [P] {X ± '} 
or 

|C<}> CX 

[S] [N] [P] « ' C<J> SA 

• S<f> 

The cosines C <f> CX etc. are available from RIGEM and [S] [N] [P] is available 
from EARRO. 

An equivalent but simpler calculation of the effect due to the Earth's 
figure is found by referring these calculations to the {y^} frame initially. 

Thus , 

(55) 

(56) 

In these equations, <J , X' are referred to {y^} and c 20 is a constant. 

Verification of RIGEM. The translational motion of the centers of mass of 
the planets and the Moon have been compared with the results given in Reference 


{ (5 sin 2 <j> - 1) C<£' CX' 
(5 sin 2 <j> - 1) C<t>' SX' • 
(5 sin 2 - 3) 

(F } = [R] T {F } . 
x y i 



A A . 

C<j> CAj 

A A ! 

C(|) SA > = 

S$ I 



[6] . Tables 10 and 11 present the results of this comparison. The comparison 
was made at J. D. 2442000.5 after an 800-day integration. Table 10 provides 
the comparison results for all planets using the 11th order option (NOR = 11) 
in RAI9s. The relatively large errors for the inner planets are due to the 
exclusion of relativity in this model. Table 11 presents comparisons for the 
Earth and Moon only for several RAI9S orders. Note that 22 is an accuracy 
parameter as discussed in Reference [7] . 

The coupled rotational-translational portion of RIGEM has been verified by 
fitting the output of RIGEM to the output of ANEAM0 as regards the physical 
libra tions in node, inclination, and longitude. Details of this comparison are 
presented in the next section. 

III. PARAMETER ESTIMATION METHOD AND PROGRAM ESTEM 

In order to verify the numerically integrated lunar librations produced by 
RIGEM and the analytic librations calculated by ANEAM0 in subroutine LURE, a 
comparison of these two approaches has been undertaken. Also, the eventual use of 
the unified model in the reduction of LURE data will necessitate a comparison of 
the observations with the model. 

The above comparisons can only be made if the proper set of initial conditions 
and model parameters is used in the comparison. The numerically integrated 
librations comprise the model in this study. For a comparison of the model with 
the analytic theory a set of initial Euler parameters and rates might be taken 
from ANEAMO since these quantities are calculated there. It has been determined 
that these parameters and rates are not accurate: enough to give the best comparison 
of the model to the libra tional theory. There are several reasons for this^ 
among them are the fact that 1) the rates are generated by an approximate numerical 
method, i.e., by differentiation of a cubic spline fit to the parameter values and 
2) the calculation of the Euler parameters themselves has approximations since the 
calculation makes use of the truncated form of the lunar theory. 

It was therefore decided to use a traditional iterative weighted least squares 
estimation of the initial state and the parameters a y 8 , Y to insure a 
"best" fit of the model to the analytical theory. 

A complication arose in this process however since the initial Euler para- 
meters and the quantities a , 8 , Y are not all independent. An iterative 


34 


Table 10. Comparison of RIGEM with Reference [6] 



|Ax| 

(Av) 

Planet 

|Ai| 

(Av/day) 


Mercury 


Venus 


Earth 


Moon 


Mars 


Jupiter 


Saturn 


Uranus 


Neptune 


Pluto 


2 x icr 5 
5 x io- 7 

1 x 10“ 10 

1 x IO -10 

4 x io' 9 

5 x icr 11 

4 x IO" 10 

2 x io -10 

2 x io -8 
2 x IO" 10 

4 x IO” 11 
1 x io" 12 

1 x IO -11 
9 x 10“ 14 

1 x 10 _1 ° 
1 x io -18 

5 x 10-1 1 
9 x 10" 14 

9 x io -11 
1 x io' 13 


4 x I0“ b 
1.3 x io” 6 

3 x io~ 9 
7 x 10" 11 

3 x io -9 
6 x IO"" 1 ! 

3 x 10“ 9 

1 x 10“ 9 

5 x io“ 8 

4 x io" 10 

2 x 10“1° 

6 x io -14 

7 x lo-H 

1 x 10~1 3 

3 x 10-H 

2 x 10 “1 3 

1 x 10“1° 

2 x 10 -I 3 

1 x 10~H 

2 X IQ-1 3 


1) RAI9S used with NOR = 11 and LL = 8. 

2) Not all forces in Reference [6] were modeled — only 
those discussed in regard to subroutines FORCE and 
FORTOR. 








Table 11. Comparison of RIGEM with Reference [6] (Continued). 


Planet 

Ax (Av) 

AY (Av) 

AZ (Av) 

Conditions 

Earth 

-1.81 x 10“ 6 

8.66 x io -7 

3.49 x IO -7 

NOR = 7 , LL = 5 

Moon 

2.75 * 10“ 5 

2.27 x lo~ 6 

3.90 x io -6 

OMPL (1) » 1 

No Earth or Lunar 
figure 

Earth 

-1.46 x 10“ 6 

8.83 x 10“ 7 

3.9 x IO” 7 

NOR = 15 , LL = 10 

Moon 

-1.02 x io~ 6 

8.92 X io' 7 

. ' 

.# 

4.9 x io “ 7 

OMPL (1) =.l 

No Earth or Lunar 
figure 

Earth 

-0.1 x 10“ 8 

-2.1 x 10“ 9 

2.2 x IO” 9 

NOR = 7 , LL - 10 

Moon 

4.7 x icr 7 

6.6 x 10~ 9 

1.7 x 10“ 7 

All Planets 

No Earth or Lunar 
figure 

Earth 

-3.87 x icr 9 

-2.69 x io -9 

-4.8 x IO -10 

NOR = 15 f LL = 10 

Moon 

1.3 x icr 8 

-1.95 x io -9 

1.44 x 10“ 9 

All Planets 

Earth figure effect 
on Lunar orbit 

Earth 

-3.72 x io -9 

-2.6 x io~ 9 

-1.4 x IO -9 

NOR = 11 , LL = 8 

Moon 

4.4 x 10“ 10 

-2.5 x 10~ 9 

-7.7 x IO -10 

All Planets 

Both Earth and Lunar 
figure effect on 
Lunar orbit 
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weighted least squares (IWLS) method was therefore programmed accounting for the 
fact that exact constraints must be satisfied between certain estimated variables. 

The Euler parameters and rates generated in ANEAMO were utilized as 
initial guesses to the IWLS process, 

A. Iterative Weighted Least Squares with Constraints 

Reference {17] provides a formulation of the iterative weighted least 
squares method when exact constraints are present, A simplified version of 
that formulation is presented here. 

The vector of observations y can be related to the vector of theoretically 

o 

calculated values from the model by 

y o = y c (X) + e (57) 

where x is a set of parameters and initial conditions and e is a vector of 
measurement errors. In this case, the y c is nonlinear in the parameters x 
and a linearization is made about a nominal set of parameters x° , viz. 


= y c (?) 


+ A (x - x°) + e 


where 


(58) 


A. . = 
11 


Sy 


c. 

x 


3x . 
3 




The traditional derivation of the best estimate, £ , to x so that 
is minimized in a weighted least squares sense is given by 


x = x° + <A T WA) - ” 1 A T W (y - y (j? 5 )) i 

O v 


(59) 


where 


T 

( ) is the transposed matrix 
( ) —1 is the inverse matrix 
W is the weighting matrix. 
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Also, if e are samples of zero-mean Gaussian independent random variables 
and if each observation is weighted with its associated error variance the 
covariance matrix of the error in x is 

P = (A T WA)' 1 , 


the standard error of the estimate of x^ is 


a. 

i 


/pt: 

XI 


and the correlation coefficient between errors in estimates of x. and x, is 

i 3 


P. . 
ID 





as shown in Reference [18]. 

The above development is modified if the parameters are not ail independent 
Reference [17] provides the following algorithm in that cases 

i) Define 

•>-> -HB 

where x is a set of "solve-for” parameters and S is a set of "exactly 
constrained' 1 parameters. 

ii) Define the exact constraints by 

f ± (x, S, BL) = 0 (i = 1, . . . , n) 

where N. is a set of n constants. • 

i 

iii) Designate one parameter from each exact constraint as a constrained 
parameter and solve for it as follows* 
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( s l ( x)\ 


is} = 


: 

IS (x) 
> n 


and determine the following partial derivatives 


[S ] - 
x 


3S i 
3xi 


3S 

n 

L 3x i 


3Sj 

3x 

n 


3S 

n 

3x 


n, 


iv) Form the residual (observed-calculated) vector 


{R} = I - - - 

I 'V 

' x - X 

where 

A 

Z are observables 
Z are computed observables 
x are a priori parameter estimates 
x are estimated values - 

v) Form the A matrix 
[A] = 




where 


1 
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-x— = A + A S 
dS x S X 


vi) Define 


J = (A + A S ) W (A + A S ) 
X S X X S X 


where W is a weighting matrix 

vii) Then, using a Newton-Raphson integration where x 
estimate of x and x^ n+ ^V is the (n + l)st estimate, there is obtained 


(n) 


is the nth 


(n+1) 

x 


(n) 


+ J 


[ (A + A 


S ) 

X 


T 


W (Z 


(X (n) )] 


(60) 


s i (n+ l) = s ± (x (n+1) ) 


(61) 


B. Program ESTEM (Version 01)- 

Program ESTEM was prepared to fit the output from RIGEM to that of ANEAMO 
by adjusting initial conditions and certain physical parameters in an iterative 
weighted least squares sense. 


Some general features of ESTEM are discussed below. 
Program ESTEM links several subprograms as follows % 
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Subroutines CTASOS and PLOT are LaRc library subroutines for inversion of 
matrices and plotting. The other subroutines have been discussed previously. 

Subroutine FORCE and FORTOR have been modified for use in ESTEM from 
the version described earlier in this report. The changes in summary form 
are : 

1. The segment that calculates forces on the Sun has been removed. 

2. An option has been included that allows the integration of only the 
Earth and Moon's translational motion (<J>MPL (1) = 2). 

3. The capability of reading the JPL ephemeris tape DE69 has been added 
for those planets not integrated. 

4. The segment that calculates the rotational motion of the Earth has 
been removed. 

5. Relativity perturbations on the planet's orbits have been included 
using the modified one-body Eddington/Roberson equations of [20} . 

6. The following forces and torques have been added: 

i) Lunar torques due to. Earth acting on all remaining 4th degree 
lunar harmonics. 

ii) Lunar torques due to interaction of C 2 2 with all second degree 
lunar harmonics. 

iii) Acceleration of lunar mass center due to Sun/Earth figure interaction 
and to Sun/Moon figure interaction. 

iv) Acceleration of Earth mass center due to Earth and lunar figures, 
due to Sun/Earth figure interaction and due to Sun/Moon figure interactions. 

The derivation for items i) and ii) follows that given in [1] while the 
derivation of items iii) and iV) follows [63 . 

This program basically performs a pre-set number of iterations based on 
equations (58) and (59) . The final estimate of the parameters is used with 
the model to compute a final set of residuals which are then plotted. 

More explicitly. 


i h 


{x} = 


00 

6i 

02 

00 

01 
02 
e 

Y 


{S> = 



( 62 ) 


IS x 3 = 


fl = 1 - 6 0 2 - 

3i 2 - 32 2 - 3 3 2 = 0 





f 2 = So 00 + 01 

3x + -32 02 + 03 03 = 0 , 

► r 



(63) 

f 3 = a (1 - 3y) 

- 3 + Y = 0 J 

> 





Si = 3 3 = ± /i~ 

- 0 O 2 - 3i 2 - S 2 2 





s 2 - S 3 - - ^ 

(00 00 + 3i 3x + 32 02^ 

» 



(64) 

Sa _ 0 .7 Y , 
3 1 - Sy 

| 






i0o/Si iSx/Si 

±3 2 /Si 0 0 

0 

0 

0 


“00/03 

“02/03 “3o/03 “3o/03 

“02/03 

0 

0 

, (65) 

0 0 

0 0 0 

0 

<1 

K 2 . 


<1 = (1 “ y 2 )/ (1 
k 2 = (3 2 1)/(1 

“ 0Y) 2 
- Sy ) 2 




(66) 


The weighting matrix is taken to be the identity matrix and the a priori 
parameter estimates and estimated values are set equal to zero in {r} „ 

The matrix of observables is 
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where p , cr , t are the physical librations in inclination, node, and 
longitude and N is the number of observations and the matrix [A] is generated 
by the "variant trajectory" method. 

r A 

In this application, the {ZJ- observables are calculated from Eckhardt’s 
theory as evaluated in ANEAMO. The calculated values {z} are calculated in 
RIGEM. 

The model used is summarized below: 
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( 69 ) 


I ! I 1 \ 

iii) IPP] = [EC] T (P] [t] T [C] T 

iv) % = tan*" 1 (-PP (3,1) /-PP (3,2)) 

6 = tan" 1 (VPP (3,1) 2 + PP (3,2) ^/PP (3,3)) ► 

ip = tan” 1 (-PP (1,3)/PP (2,3)) 


( 70 ) 


v) n = a a + b a fc + V (71) 

<r = a (f + + c^t 2 

vi) p = 0 - I 

a * ip - S2 (72) 

T = - n - (f + $ 


In the above, equations (65) and (66) are derived in Reference [1], 
equations (67) and (68) were discussed earlier as equations (43) and (44). 
Equations (69) are available in Reference [2] , and equation (70) was discussed 
earlier as equation (45) . 

C. A Verification Run 

The results of the verification run are shown in Figure 3. There the 
residuals in p , la , and r in arc seconds are shown versus time in days 
past the epoch of J. D. 2441200.5. The fit was made over 1100 days. 

The maximum residuals in p and la are about ±3" and those in T are 
about ±10" . This plot indicates that no gross errors exist but some subtle 
inconsistencies still exist in the comparison. The residuals in all angles 
should be periodic with maximum amplitudes of <1’.'5 . Work is continuing on 
improvement of these residuals as well as a comparison with a more extensive 
version of the theory (Eckhardt's 300 series). 

The following conditions applied for the comparison: 

ESTEM: 

i) Initial conditions for planetary motions per Reference .['61. 

ii) Planetary masses per Reference [6] . 

iii) Lunar nominal Euler parameters and rates 









I 


-.9805662538581 

.1923761339498 

3.629025185335 -2 


So 
01 
02 

g 3 = -1.281649334705 


-2 


iv) a = 4.023"** 
3 = 6.293“ 4 


v = 2,27 


-4 


00 

01 

02 

03 


1. 872186798183*7 3 
7. 693100373778 -3 
1.478710744051" 2 
1. 410632074475 ~ 2 


lunar inertia ratios. 


v) Sampling interval = 3 days, 

vi) Integration order - H (ll = 8) in RAI9S. 

vii) Only second degree lunar harmonics considered. 

viii) C 20 = -1-082637* 3 Earth Oblateness. 




ix) Speed of light = 299792.5 km/sec. 

x) Relativity perturbations calculated using Eddington/Robertson form. 
Reference [20] . 

xi) No tidal coupling effect on Earth and lunar orbits, 

xii) Earth and Moon integrated with remaining planetary motions read 
from JPL DE69 ephemeris tape. 

ANEAMO: 

i) Eckhardt j = 3 model with coefficients as listed in Appendix A. 

IV. FUTURE WORK 

The next phases of work under NSG-1152 will involve the following: 

1. Development of observational equations (laser ranging data) and partial 
derivatives ; 

2. Adopting a solution parameter set to be estimated from data; 

3. General optimization in terms of running time and storage requirements 
of all programs; and 
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4. Development of variational equations for use in parameter estimation 
process. 

To implement the above a series of computer programs based on ANEAMO, 

RIGEM, and ESTEM is envisional. Collectively these programs are referred to 
as EMSYS. 

The individual programs will have the following capabilities: 

EMDYN. i) Numerically integrate all translational, rotational and 
variational equations of motion from a specified set of initial conditions and 
parameters; ii) Calculate numerical partials of some quantities if necessary; 
and iii) Output "ephemeris" of all quantities on tape 10. 

EMOBS. i) Read LURE data tape; ii) Apply any necessary corrections to 
the data; and iii) Write resulting observations on tape 11. 

EMNORM. i) Reads tapes 10 and 11; ii) Interpolates tape 10 at observation 
times; iii) Computes relativistic time delay; iv) Calculates all partial 
derivatives using analytical results and data from the variational equations 
as supplied on tape 10; v) Forms the normal equations; and vi) Outputs 
residuals and normal equations on tape 12. 

EMEST. i) Reads tape 12; ii) Solves normal equations for best estimate 
of parameters; iii) Outputs useful information in printed and plotted form; 
and iv) Computes covariance matrix and other statistical quantities. 
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APPENDIX A 


ECKHARDT'S LIBRATION THEORY 
(J = 3 ) 
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APPENDIX B 


ERRATA FOR FIRST 
SEMIANNUAL PROGRESS REPORT 
ON NSG-1152 


* 

Corrected Pages 
from Reference [1] Follow 



JS 


Inverting these expressions provides 


r = A\ + A* + a| 


X - tan ** 1 (A2/A1) , 


and 


( 41 ) 


<p = tan ” 1 (A 


\/M + 4 ) 


Since the unit vectors k. are related to those of the spherical 
polar system by 


ki 

Jc 2 


= -X. 


= -X, 


( 42 ) 


U = 


the inertial angular velocity of the axes {Z^} can be written as 


9 • 


W Z/X' - A + ( I )= ^ I 3 ~ < PX^ 


( 43 ) 


The abdve vector may be projected on the {Z^} axes providing 


“z/X' = * [ cos ^ ^3 “ sin <J> &i] + 5> *2 


( 44 ) 


The components 


Wj = 
* 

t *?2 = 

i 

W3 = 


-X sin <j> 

- 9 ■ 

♦ 

X cos <j> 


treated together. Finally the term U3 and the remaining terms in 
U4® ^ will be treated [see eqs. (86) and ( 87 )]. 

The reference axes for the second order and coupling terms 
are {y^} . Thus the c!^ and s!^ are functions of the orien- 
tation angles. These functions are 


C20 = 


_1 f T yiy\ f x y2yz . z S 

.■ 2 M' L 2 y 3 y 3J 


1 r a + b + 

2 m' L 2 


- - | (A'y 2 + B'y' 2 + C'Y 


..2)j 


= 17^*3 


a ' 2 M' 


[ctyA' + a ' y 'B’ + o' ' y* 'C' ] 


e t = 1- T * 

521 a .2 M . W 


( 80 ) 


a' 2 M* 


- Ly 3 A' + y' 3 'B' + Y* ’ 3' 'c'3 


C22 “ 


4 a 


- [A' 

' 2 m* - 


' (B 2 - a 2 ) + B' ( 3 ' 2 - a' 2 ) 


+ C' (8' ’ 2 


a 1 ' 2 )j 


S 2 2 = 


4 a' 2 M* 


[a $A' + a 1 S'B' + a' ' 3' 'C* ] 




M = GMr“ 5 a 2 (B' - C\) |6P 40 jc 20 

Z J, 1 ' 


^6P 40 ) c 20 (a 1 *y* + y* ' a' ) 

+ c^fB'S'* - a'a ,, )}+ 3P 4 1 { -c 2 2 ( ot ' 1 y ' + y' 'a')cX 
- sX(y' 3'' + Y * * 3 ' ) (c 2 0 + c 22 )( 

3 ' 3 ' ' ) + s2X (a ' ’ 3 ' + a’3’ ’)} 


P 42 . 

— ^ — c 20 | c2X (a'a 1 ' 

P 43 

3P 4 2 c 22 c 2Xy ' Y ’ 1 - — y 


«■/*• - __ c 22 |e3X(a ,, Y' + a'Y'') 

p 4 4 


* H*f 

+ s3X<3"y' + B'y”)} - -y c 22 {c4X{3 , 3’' “ a'a") 

+ s4X(a"3' + 3' 'a')} J 

M 22 = GMr “ 5 a 2 (C* -A') £ 6P 4 0 { c 2 0 (y ' ’ a + yet " ) 

+ c 22 (33 ’ ' - aa* ' ) } + 3P 41 { -c 22 (ay' 3 + ya ' ' )cX 
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